O bloco a seguir carrega as dependências necessárias. Estão são:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker
import IPython
import scipy.signal
import scipy.stats
from scipy.fftpack import fft, ifft, rfft
from scipy.io.wavfile import read as audioread
from scipy.io.wavfile import write as audiowrite
import seaborn as sns
E abaixo a função de geração do TSP do Hani:
from gen_tsp2 import gen_tsp2
Uso:
tsp, t = gen_tsp2(T=1200, fs=44.1, bw=22.05, bs=0, ta=120, tb=50, plot=False)
Parâmetros:
| Parâmetro | Descrição | Unidade | Valor Padrão |
|---|---|---|---|
| T | Duração do TSP | ms | 1200 |
| fs | Taxa de amostragem | kHz | 44.1 |
| bw | Largura da faixa de passagem | kHz | 22.05 |
| bs | Deslocamento da faixa de passagem | kHz | 0 |
| ta | Início do TSP | s | 120 |
| tb | Taxa de crescimento do atraso de grupo | ms/kHz | 50 |
| plot | Mostrar os gráficos do TSP obtido | - | False |
Valores de retorno:
| Valor | Descrição |
|---|---|
| tsp | TSP no domínio do tempo |
| t | Vetor temporal associado a x |
O bloco a seguir gera um TSP de 0 a 22,05 kHz com duração de 5s e o salva como um arquivo .wav:
# Parâmetros do TSP
T = 5000
fs = 44.1
bw = fs / 2
bs = 0
ta = T / 10
tb = 200
# Atribui o TSP gerado à variável 'tsp' e o vetor temporal associado a 't'
# Troque para plot=True se quiser ver o gráfico do TSP
tsp, t = gen_tsp2(T, fs, bw, bs, ta, tb, plot=False)
t /= 1000
# Salva o TSP gerado como um arquivo .wav
audiowrite(filename='tsp_python.wav',
rate=int(fs*1e3),
data=tsp.astype(np.float32))
Após reproduzir o TSP e gravar a resposta da sala, os arquivos produzidos são carregados:
# A lista abaixo contém os nomes dos arquivos a serem carregados, entre
# aspas e separados por vírgulas
lista_arquivos = ['exportadas/P1Mic1g80 - 2017-06-03 - Nossa Senhora da Conceição - 7s.wav',
'exportadas/P1Mic2g80 - 2017-06-03 - Nossa Senhora da Conceição - 7s.wav',
#'exportadas/P1Mic3g80 - 2017-06-03 - Nossa Senhora da Conceição - 7s.wav',
#'exportadas/P1Mic4g80 - 2017-06-03 - Nossa Senhora da Conceição - 7s.wav',
#'exportadas/P1Mic5g80 - 2017-06-03 - Nossa Senhora da Conceição - 7s.wav',
#'exportadas/P1Mic6g80 - 2017-06-03 - Nossa Senhora da Conceição - 7s.wav',
]
# Após o bloco "for" a seguir, a variável 'respostas' conterá uma lista
# das respostas ao tsp carregadas
respostas = []
for arquivo in lista_arquivos:
rate, resp = audioread(arquivo) # carrega a resposta do arquivo .wav
resp = resp / np.max(np.abs(resp)) * .98 # normaliza para estar entre -1 e 1
respostas.append(resp[:len(t)]) # põe na lista de respostas
# num_respostas é o número de respostas carregadas
num_respostas = len(respostas)
# players de áudio para ouvir as respostas carregadas
audio0 = IPython.display.Audio(data=tsp, rate=(fs * 1e3)) # tsp
label = IPython.display.Markdown('TSP')
IPython.display.display(label, audio0)
for indice, resposta in enumerate(respostas):
titulo = IPython.display.Markdown(data='Resposta {}'.format(indice + 1)) # título
audio = IPython.display.Audio(data=resposta, rate=(fs * 1e3)) # áudio
IPython.display.display(titulo)
IPython.display.display(audio)
As respostas obtidas no domínio do tempo estão ilustradas abaixo:
# cria gráficos suficientes para plotar o TSP e todas as respostas
fig, graficos = plt.subplots(num_respostas + 1, figsize=(15, 20)) # 'graficos' é uma lista de gráficos
plt.tight_layout() # apenas estético
# plota o TSP no primeiro gráfico (0)
graficos[0].plot(t, tsp) # plota o tempo no eixo horizontal e tsp no vertical
graficos[0].set_title('TSP') # título do gráfico do TSP
# plota cada resposta nos graficos seguintes
for indice, resposta in enumerate(respostas): # para cada resposta
graficos[indice + 1].plot(t, resposta) # plota o gráfico da resposta
graficos[indice + 1].set_title('Resposta {}'.format(indice + 1)) # título da resposta
# insere o texto dos eixos
fig.text(.5, .001, 'Time (ms)', ha='center')
fig.text(-.02, .5, 'Normalized Amplitude', va='center', rotation='vertical')
#fig.savefig('Respostas_no_tempo.pdf', format='pdf')
Abaixo é feita a FFT de cada resposta:
# calcula um vetor de frequências 'f' (análogo ao vetor 't')
f = np.arange(0, len(t), 1) * (fs*1000) / len(t)
# calcula a FFT das respostas e as coloca em uma lista
fft_respostas = [fft(resp) for resp in respostas]
# e as respostas em frequência resultantes são passadas pra decibel (20 x log10 (abs(x))
db_fft_respostas = [20 * np.log10(np.abs(resp)) for resp in fft_respostas]
As respostas em frequência obtidas estão plotadas abaixo:
# cria uma lista das bandas padrão de 1/3 de oitava em kHz
freqs = np.array([.020, .025, .032, .040, .050, .063, .080, .100, .125, .160, .200, .250, .315, .400, .500, .630, .800,
1.0, 1.25, 1.6, 2.0, 2.5, 3.15, 4.0, 5.0, 6.3, 8.0, 10.0, 12.5, 16.0])
# cria gráficos suficientes para plotar o TSP e todas as respostas
fig, graficos = plt.subplots(num_respostas + 1, figsize=(15, 20)) # 'graficos' é uma lista de gráficos
plt.tight_layout() # apenas estético
# calcula a resposta em frequência do TSP para plotar (deve ser uma reta)
fft_tsp = fft(tsp)
db_fft_tsp = 20 * np.log10(np.abs(fft_tsp))
# plota o TSP no primeiro gráfico (0)
graficos[0].plot(f, db_fft_tsp) # plota as frequências no eixo horizontal e o tsp no vertical
graficos[0].set_title('TSP') # título do gráfico do TSP
graficos[0].set_xscale('log', basex=2)
graficos[0].set_xticks(1000 * freqs)
graficos[0].xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
graficos[0].axis([20, (fs*1000)/2, 0, 60]) # ajusta os limites do gráfico
# plota cada resposta nos graficos seguintes
for indice, resposta in enumerate(db_fft_respostas): # para cada resposta em frequência
graficos[indice + 1].plot(f, resposta) # plota o gráfico da resposta
graficos[indice + 1].set_title('Resposta {}'.format(indice + 1)) # título da resposta
graficos[indice + 1].axis([20, (fs*1000)/2, 0, 60]) # ajusta os limites dos eixos
# escala logarítmica
graficos[indice + 1].set_xscale('log', basex=2)
graficos[indice + 1].set_xticks(1000 * freqs)
graficos[indice + 1].xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
# põe linhas em todos os dós e fás
# for i in range(0, 10):
# c_freq = 32.7032 * 2**i
# f_freq = 43.6535 * 2**i
# graficos[indice + 1].axvline(c_freq, color='red', ls='-', label='Dó')
# graficos[indice + 1].axvline(f_freq, color='black', ls='-', label='Fá')
# insere o texto dos eixos
fig.text(.5, .001, 'Frequency (Hz)', ha='center')
fig.text(-.02, .5, 'Amplitude (dB)', va='center', rotation='vertical')
#fig.savefig('Resposta_em_frequencia.pdf', format='pdf')
A função 'get_ir_from_fr' abaixo calcula a resposta ao impulso da sala a partir de uma resposta em frequência. O cálculo é feito no domínio da frequência, i.e., se a resposta da sala que se deseja obter é $H$, o TSP é $X$ e a resposta gravada é $Y$, pode-se obter $H$ por $$H=\frac{Y}{X}=YX^{-1}$$
def get_ir_from_fft(fft_resp, fft_tsp):
fft_ir = fft_resp / fft_tsp # deconvolui o TSP do sinal gravado
ir = np.real(ifft(fft_ir)) # obtém a resposta ao impulso através da transformada inversa
ir = ir / np.max(np.abs(ir)) * 0.98 # normalização
return ir
Uso:
ir = get_ir_from_fft(fft_resp, fft_tsp)
Parâmetros:
| Parâmetro | Descrição |
|---|---|
| fft_resp | FFT da resposta da sala ao TSP |
| fft_tsp | FFT do TSP utilizado |
Valores de retorno:
| Valor | Descrição |
|---|---|
| ir | Vetor com a resposta ao impulso calculada |
A linha a seguir calcula as respostas ao impulso a partir das respostas em frequência:
irs = [get_ir_from_fft(fft_resp, fft_tsp) for fft_resp in fft_respostas]
As respostas ao impulso obtidas estão plotadas abaixo:
# cria gráficos suficientes para plotar o TSP e todas as respostas
fig, graficos = plt.subplots(num_respostas + 1, figsize=(15, 20)) # 'graficos' é uma lista de gráficos
plt.tight_layout() # apenas estético
# plota o TSP no primeiro gráfico (0)
graficos[0].plot(t, tsp) # plota as frequências no eixo horizontal e o tsp no vertical
graficos[0].set_title('TSP') # título do gráfico do TSP
# plota cada resposta nos graficos seguintes
for indice, resposta in enumerate(irs): # para cada resposta ao impulso
graficos[indice + 1].plot(t, resposta) # plota o gráfico da resposta
graficos[indice + 1].set_title('Resposta {}'.format(indice + 1)) # título da resposta
graficos[indice + 1].axis([0, 6, -1, 1])
# insere o texto dos eixos
fig.text(.5, .01, 'Time (ms)', ha='center')
fig.text(-.02, .5, 'Normalized Amplitude', va='center', rotation='vertical')
#fig.savefig('Resposta_ao_impulso.pdf', format='pdf')
# players de áudio para ouvir as respostas calculadas
for indice, resposta in enumerate(irs):
titulo = IPython.display.Markdown(data='Resposta {}'.format(indice + 1)) # título
audio = IPython.display.Audio(data=resposta, rate=(fs * 1e3)) # áudio
IPython.display.display(titulo)
IPython.display.display(audio)
A partir das respostas ao impulso obtidas, é possível calcular o tempo de reverberação para cada faixa de frequências. O cálculo é feito da seguinte forma:
def rt60(ir, freqs, fs, beg=-5, end=-25, mul=3, plot=False):
rt60 = [] # lista que conterá o RT60 para cada frequência pedida
if plot:
fig, axs = plt.subplots(len(freqs))
if len(freqs) == 1:
axs = [axs]
for index, freq in enumerate(freqs): # para cada frequência
# calcula as frequências de corte do filtro
freq *= 1000 # frequência em Hz
low = freq / 2**(1/6) # 2**(1/6) -> 1 terço de oitava
low = low / (fs/2 * 1e3) # normalização
high = freq * 2**(1/6)
high = high / (fs/2 * 1e3)
# calcula os parâmetros do filtro
num, den = scipy.signal.butter(N=3, Wn=[low, high], btype='bandpass') # passa-faixas de terceira ordem
# filtra o sinal e o normaliza
filtered = scipy.signal.lfilter(num, den, ir)
filtered = np.abs(filtered) / np.max(np.abs(filtered))
# obtém a representação do sinal em dB para encontrar o limite de integração
db_signal = 20 * np.log10(filtered)
# obtém o início do impulso
impulse = np.where(filtered == np.max(filtered))[0][0]
# encontra o ponto em que o sinal atravessa o nível de ruído
window = int(fs * 50) # 50 ms
avg_signal = np.convolve(db_signal, np.ones(window) / window, mode='valid') # média móvel
# estende o vetor da média móvel
padding = int((len(db_signal) - len(avg_signal)) / 2)
avg_signal = np.pad(avg_signal, padding, 'edge')
threshold = int(len(db_signal) * 0.8) # último décimo do sinal
noise_level = np.average(avg_signal[threshold:]) + 5 # calcula o nível de ruído como a média dos 10% finais
t1 = np.abs(avg_signal[impulse:] - noise_level).argmin() + impulse # cruzamento do nível de ruído
# realiza a integração de Schroeder
sch = np.cumsum(filtered[t1:0:-1] ** 2)[::-1] # sinal integrado
sch_db = 10 * np.log10(sch / np.max(sch)) # sinal integrado em dB
# Encontra os pontos de cruzamento para a regressão
#beg = -5 # inicio em -5 dB (ISO 3382-1)
#end = -25 # -25 dB para RT20
#mul = 3 # 3*RT20 = RT60
first = sch_db[np.abs(sch_db - beg).argmin()] # cruzamento com -5 dB
first_sample = np.where(sch_db == first)[0][0] # sample em que o cruzamento ocorre
last = sch_db[np.abs(sch_db - end).argmin()] # cruzamento com -25 dB
last_sample = np.where(sch_db == last)[0][0]
# realiza a regressão linear nos pontos encontrados
time = np.arange(first_sample, last_sample + 1) / fs
sig = sch_db[first_sample:last_sample + 1]
slope, intercept, r_value, p_value, std_err = scipy.stats.linregress(time, sig)
reg_init = (beg - intercept) / slope # sample em que a regressão cruza -5 dB
reg_end = (end - intercept) / slope # sample em que a regressão cruza -25 dB
rt60.append(mul * (reg_end - reg_init))
if plot:
line = slope*np.arange(0, last_sample / fs + 1000) + intercept
t = np.arange(0, len(db_signal)) / fs
axs[index].plot(line)
axs[index].plot(t, db_signal)
axs[index].axvline(x=reg_init)
axs[index].axvline(x=reg_end)
axs[index].axhline(y=beg)
axs[index].axhline(y=end)
plt.show()
return rt60
Uso:
lista_rt60 = rt60(ir, freqs, fs, plot=False)
Parâmetros:
| Parâmetro | Descrição | Unidade | Valor Padrão |
|---|---|---|---|
| ir | Resposta ao impulso da sala | - | - |
| freqs | Lista com as frequências desejadas | kHz | - |
| fs | Taxa de amostragem | kHz | - |
| plot | Exibir gráfico para cada frequência | - | False |
Valores de retorno:
| Valor | Descrição | Unidade |
|---|---|---|
| lista_rt60 | Lista com o RT60 para cada frequência | ms |
O bloco a seguir calcula o RT60 as bandas padrão de 1/3 de oitava e os exibe em um gráfico, para cada resposta ao impulso obtida.
# cria uma lista das bandas padrão de 1/3 de oitava em kHz
freqs = np.array([.063, .080, .100, .125, .160, .200, .250, .315, .400, .500, .630, .800,
1.0, 1.25, 1.6, 2.0, 2.5, 3.15, 4.0, 5.0, 6.3, 8.0, 10.0, 12.5, 16.0])
# calcula e plota os RT60 para cada frequência
lista_edt = list()
lista_rt60_por_rt20 = list() # rt60 calculado para cada amostra
lista_rt60_por_rt30 = list()
lista_rt60_por_rt60 = list()
for resposta in irs: # para cada resposta ao impulso
lista_edt.append(np.array(rt60(resposta, freqs, fs, beg=0, end=-10, mul=1))) # obtém o edt
lista_rt60_por_rt20.append(np.array(rt60(resposta, freqs, fs, beg=-5, end=-25, mul=3))) # obtém o RT60 pelo RT20
lista_rt60_por_rt30.append(np.array(rt60(resposta, freqs, fs, beg=-5, end=-35, mul=2))) # obtém o RT60 pelo RT30
lista_rt60_por_rt60.append(np.array(rt60(resposta, freqs, fs, beg=-5, end=-65, mul=1))) # obtém o RT60 diretamente
# cria uma lista das bandas padrão de 1/3 de oitava em kHz
freqs = np.array([.063, .080, .100, .125, .160, .200, .250, .315, .400, .500, .630, .800,
1.0, 1.25, 1.6, 2.0, 2.5, 3.15, 4.0, 5.0, 6.3, 8.0, 10.0, 12.5, 16.0])
# cria um número suficiente de gráficos
fig, axs = plt.subplots(num_respostas + 1, figsize=(15, 20)) # um gráfico adicional para a média
plt.tight_layout()
if num_respostas == 1:
axs = [axs]
# calcula e plota os RT60 para cada frequência
for indice in range(len(irs)): # para cada resposta ao impulso
axs[indice].plot(freqs*1000, np.array([
#lista_edt[indice] / 1000,
lista_rt60_por_rt20[indice] / 1000,
#lista_rt60_por_rt30[indice] / 1000,
#lista_rt60_por_rt60[indice] / 1000,
]).transpose()) # plota os tempos obtidos em um gráfico
axs[indice].set_title('Resposta {}'.format(indice + 1))
axs[indice].axis([64, fs*1000/2, 0, 7]) # alterar os limites de acordo com a situação
# linha mediana
mediana = np.median(lista_rt60_por_rt20[indice])
rt60_1k = lista_rt60_por_rt20[indice][12] /1000
axs[indice].axhline(y=rt60_1k, color='red', ls='--', label='RT60 em 1.000 Hz ={:.2f} s'.format(rt60_1k))
# escala logarítmica
axs[indice].set_xscale('log', basex=2)
axs[indice].set_xticks(1000 * freqs)
axs[indice].xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
# legenda
handles, labels = axs[indice].get_legend_handles_labels()
axs[indice].legend(handles, labels)
# calcula a média e o desvio padrão das amostras em cada frequência
media_rt60, desvio_rt60 = zip(*[(freq.mean() / 1000, freq.std() / 1000) for freq in np.array(lista_rt60_por_rt20).transpose()])
#axs[-1].plot(freqs*1000, np.array(media_rt60))
axs[-1].errorbar(freqs*1000, media_rt60, yerr=desvio_rt60, ecolor='black', label='RT60 com desvio padrão')
axs[-1].set_title('RT60 médio da igreja')
axs[-1].axis([64, fs*1000/2, 0, 2.5]) # alterar os limites de acordo com a situação
# linha mediana
mediana = np.median(media_rt60) / 1000
rt60_1k = media_rt60[12]
axs[-1].axhline(y=rt60_1k, color='red', ls='--', label='RT60 em 1.000 Hz = {:.2f} s'.format(rt60_1k))
# escala logarítmica
axs[-1].set_xscale('log', basex=2)
axs[-1].set_xticks(1000 * freqs)
axs[-1].xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
# legenda
handles, labels = axs[-1].get_legend_handles_labels()
axs[-1].legend(handles, labels)
# insere o texto dos eixos
fig.text(.5, .001, 'Frequency (Hz)', ha='center')
fig.text(-.02, .5, 'RT60 (s)', va='center', rotation='vertical')
#fig.savefig('RT60.pdf', format='pdf')
Os parâmetros de claridade são a razão entre a energia sonora antes e depois de um determinado instante, i.e. $$C_t=10 \log_{10}\left( \frac{ \int_0^t p^2(\tau) d\tau }{ \int_t^\infty p^2(\tau) d\tau } \right).$$ Os valores mais utilizados são 50 ms (correlacionado com a inteligibilidade da fala) e 80 ms (para música). A função abaixo calcula a claridade para qualquer valor de tempo e determinadas frequências, a partir da resposta ao impulso:
def clarity(ir, time, freqs, fs):
init_sample = np.where(ir == np.max(ir))[0][0] # sample do início do impulso
c = [] # lista onde serão armazenados os valores de claridade
for freq in freqs: # para cada frequência
# obtém as frequências de corte do filtro
freq *= 1e3 # frequencia em Hz
low = freq / 2**(1/6) # 1/3 de oitava
low = low / (fs/2 * 1e3) # normalização
high = freq * 2**(1/6)
high = high / (fs/2 * 1e3)
# calcula os parâmetros do filtro
num, den = scipy.signal.butter(3, [low, high], 'bandpass')
# filtra o sinal na faixa desejada
filtered_signal = scipy.signal.lfilter(num, den, ir)
# calcula a claridade pela fórmula acima
h2 = filtered_signal ** 2 # quadrado do sinal
t = init_sample + int(time*fs + 1) # valor de t na fórmula
c.append(10 * np.log10(np.sum(h2[init_sample:t]) / np.sum(h2[t:]))) # integração
return c
Uso:
lista_claridades = clarity(ir, time, freqs, fs)
Parâmetros:
| Parâmetro | Descrição | Unidade |
|---|---|---|
| ir | Resposta ao impulso calculada | - |
| time | Instante desejado, ex. 50 para C50 | ms |
| freqs | Lista de frequências | kHz |
| fs | Taxa de amostragem | kHz |
Valor retornado:
| Valor | Descrição |
|---|---|
| lista_claridades | Lista com o valor da claridade para o tempo escolhido em cada frequência |
O bloco abaixo calcula os valores de C50 para várias frequências das respostas obtidas e os plota em gráficos.
# cria uma lista das bandas padrão de 1/3 de oitava em kHz
freqs = np.array([.063, .080, .100, .125, .160, .200, .250, .315, .400, .500, .630, .800,
1.0, 1.25, 1.6, 2.0, 2.5, 3.15, 4.0, 5.0, 6.3, 8.0, 10.0, 12.5, 16.0])
# cria um número suficiente de gráficos
fig, axs = plt.subplots(num_respostas + 1, figsize=(15, 20))
plt.tight_layout()
if num_respostas == 1:
axs = [axs]
# plota os C50 para cada frequência
lista_c50 = list()
for indice, resposta in enumerate(irs): # para cada resposta ao impulso
lista_c50.append(clarity(resposta, 50, freqs, fs)) # obtém o C50
axs[indice].plot(freqs*1000, lista_c50[indice]) # plota os tempos obtidos em um gráfico
axs[indice].set_title('Resposta {}'.format(indice + 1))
axs[indice].axis([min(freqs)*1000, fs*1000/2, -10, 15]) # alterar os limites de acordo com a situação
# linha mediana
mediana = np.median(lista_c50[indice])
c50_1k = lista_c50[indice][12]
axs[indice].axhline(y=c50_1k, color='red', ls='--', label='C50={:.2f} dB'.format(c50_1k))
# escala logarítmica
axs[indice].set_xscale('log', basex=2)
axs[indice].set_xticks(1000 * freqs)
axs[indice].xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
# legenda
handles, labels = axs[indice].get_legend_handles_labels()
axs[indice].legend(handles, labels)
# calcula a média e o desvio padrão das amostras em cada frequência
media_c50, desvio_c50 = zip(*[(freq.mean(), freq.std()) for freq in np.array(lista_c50).transpose()])
#axs[-1].plot(freqs*1000, media_c50)
axs[-1].errorbar(freqs*1000, media_c50, yerr=desvio_c50)
axs[-1].set_title('Média')
# escala logarítimica
axs[-1].set_xscale('log', basex=2)
axs[-1].set_xticks(1000 * freqs)
axs[-1].xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
# linha mediana
mediana = np.median(media_c50)
c50_1k = media_c50[12]
axs[-1].axhline(y=c50_1k, color='red', ls='--', label='C50={:.2f} dB'.format(c50_1k))
axs[-1].axis([min(freqs) * 1000, fs*1000/2, -10, 10]) # limites
# legenda
handles, labels = axs[-1].get_legend_handles_labels()
axs[-1].legend(handles, labels)
# insere o texto dos eixos
fig.text(.5, .001, 'Frequency (Hz)', ha='center')
fig.text(-.02, .5, 'C50 (dB)', va='center', rotation='vertical')
#fig.savefig('C50.pdf', format='pdf')
E o bloco abaixo faz o mesmo com o C80 de cada resposta
# cria uma lista das bandas padrão de 1/3 de oitava em kHz
freqs = np.array([.063, .080, .100, .125, .160, .200, .250, .315, .400, .500, .630, .800,
1.0, 1.25, 1.6, 2.0, 2.5, 3.15, 4.0, 5.0, 6.3, 8.0, 10.0, 12.5, 16.0])
# cria um número suficiente de gráficos
fig, axs = plt.subplots(num_respostas + 1, figsize=(15, 20))
plt.tight_layout()
if num_respostas == 1:
axs = [axs]
# plota os C80 para cada frequência
lista_c80 = list()
for indice, resposta in enumerate(irs): # para cada resposta ao impulso
lista_c80.append(clarity(resposta, 80, freqs, fs)) # obtém o C80
axs[indice].plot(freqs*1000, lista_c80[indice]) # plota os tempos obtidos em um gráfico
axs[indice].set_title('Resposta {}'.format(indice + 1))
axs[indice].axis([min(freqs)*1000, fs*1000/2, -17.5, 15]) # alterar os limites de acordo com a situação
axs[indice].fill_between(freqs*1000, -2, 2, color='green', alpha=.25,label='Faixa de variacão ótima (-2dB à +2dB)')
axs[indice].fill_between(freqs*1000, -5, -2, color='green', alpha=.1, label='Faixa de variacão típica (-5dB à +5dB)')
axs[indice].fill_between(freqs*1000, 2, 5, color='green', alpha=.1)
# linha mediana
mediana = np.median(lista_c80[indice])
c80_1k = lista_c80[indice][12]
axs[indice].axhline(y=c80_1k, color='red', ls='--', label='C80 em 1.000Hz = {:.2f} dB'.format(c80_1k))
# escala logarítmica
axs[indice].set_xscale('log', basex=2)
axs[indice].set_xticks(1000 * freqs)
axs[indice].xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
# legenda
handles, labels = axs[indice].get_legend_handles_labels()
axs[indice].legend(handles, labels, loc='lower right')
# calcula a média e o desvio padrão das amostras em cada frequência
media_c80, desvio_c80 = zip(*[(freq.mean(), freq.std()) for freq in np.array(lista_c80).transpose()])
#axs[-1].plot(freqs*1000, media_c80)
axs[-1].errorbar(freqs*1000, media_c80, yerr=desvio_c80, label='C80 com desvio padrão')
axs[-1].set_title('C80 médio da igreja')
axs[-1].fill_between(freqs*1000, -2, 2, color='green', alpha=.25, label='Faixa de variacão ótima (-2dB à +2dB)')
axs[-1].fill_between(freqs*1000, -5, -2, color='green', alpha=.1, label='Faixa de variacão típica (-5dB à +5dB)')
axs[-1].fill_between(freqs*1000, 2, 5, color='green', alpha=.1)
# escala logarítimica
axs[-1].set_xscale('log', basex=2)
axs[-1].set_xticks(1000 * freqs)
axs[-1].xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
# linha mediana
mediana = np.median(media_c80)
c80_1k = media_c80[12]
axs[-1].axhline(y=c80_1k, color='red', ls='--', label='C80 em 1.000 Hz = {:.2f} dB'.format(c80_1k))
axs[-1].axis([min(freqs) * 1000, fs*1000/2, -17.5, 15]) # limites
# legenda
handles, labels = axs[-1].get_legend_handles_labels()
axs[-1].legend(handles, labels, loc='lower right')
# insere o texto dos eixos
fig.text(.5, .001, 'Frequency (Hz)', ha='center')
fig.text(-.02, .5, 'C80 (dB)', va='center', rotation='vertical')
#fig.savefig('C80.pdf', format='pdf')
A definição é um outro indicador da inteligibilidade da fala. É a razão entre a energia sonora nos primeiros 50 ms e a de todo o período de captura, i.e. $$D_{50} = \frac{ \int_0^{50} p^2(\tau) d\tau }{ \int_0^\infty p^2(\tau) d\tau }.$$ A função abaixo calcula a claridade de uma resposta ao impulso em várias frequências.
def definition(ir, freqs, fs):
init_sample = np.where(ir == np.max(ir))[0][0] # amostra do início do impulso
d = [] # lista onde serão armazenados os valores de definição
for freq in freqs: # para cada frequência
# obtém as frequências de corte do filtro
freq *= 1e3 # frequencia em Hz
low = freq / 2**(1/6) # 1/3 de oitava
low = low / (fs/2 * 1e3) # normalização
high = freq * 2**(1/6)
high = high / (fs/2 * 1e3)
# calcula os parâmetros do filtro
num, den = scipy.signal.butter(3, [low, high], 'bandpass')
# filtra o sinal na faixa desejada
filtered_signal = scipy.signal.lfilter(num, den, ir)
# calcula o D50 pela fórmula acima
h2 = filtered_signal ** 2 # quadrado do sinal
t = init_sample + int(50*fs + 1) # 50 ms após o início do impulso
d.append(np.sum(h2[init_sample:t]) / np.sum(h2[init_sample:])) # integração
return d
Uso:
lista_definicoes = definition(ir, freqs, fs)
Parâmetros:
| Parâmetro | Descrição | Unidade |
|---|---|---|
| ir | Resposta ao impulso calculada | - |
| freqs | Lista de frequências | kHz |
| fs | Taxa de amostragem | kHz |
Valor retornado:
| Valor | Descrição |
|---|---|
| lista_definitions | Lista com o valor da definição da resposta ir em cada frequência |
O bloco a seguir calcula a definição em várias frequências de cada uma das respostas obtidas e plota o resultado em um gráfico.
# cria uma lista das bandas padrão de 1/3 de oitava em kHz
freqs = np.array([.063, .080, .100, .125, .160, .200, .250, .315, .400, .500, .630, .800,
1.0, 1.25, 1.6, 2.0, 2.5, 3.15, 4.0, 5.0, 6.3, 8.0, 10.0, 12.5, 16.0])
# cria um número suficiente de gráficos
fig, axs = plt.subplots(num_respostas + 1, figsize=(15, 20))
plt.tight_layout()
if num_respostas == 1:
axs = [axs]
# plota os D50 para cada frequência
lista_d50 = list()
for indice, resposta in enumerate(irs): # para cada resposta ao impulso
lista_d50.append(definition(resposta, freqs, fs)) # obtém o D50
axs[indice].plot(freqs*1000, lista_d50[indice]) # plota os tempos obtidos em um gráfico
axs[indice].set_title('Resposta {}'.format(indice + 1))
axs[indice].axis([min(freqs)*1000, fs*1000/2, 0, 1]) # alterar os limites de acordo com a situação
axs[indice].fill_between(freqs*1000, .3, .7, color='green', alpha=.2, label='Faixa de variação típica (0.3 à 0.7)')
# linha mediana
mediana = np.median(lista_d50[indice])
d50_1k = lista_d50[indice][12]
axs[indice].axhline(y=d50_1k, color='red', ls='--', label='D50={:.2f}'.format(d50_1k))
# escala logarítmica
axs[indice].set_xscale('log', basex=2)
axs[indice].set_xticks(1000 * freqs)
axs[indice].xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
# legenda
handles, labels = axs[indice].get_legend_handles_labels()
axs[indice].legend(handles, labels, loc='lower right')
# calcula a média e o desvio padrão das amostras em cada frequência
media_d50, desvio_d50 = zip(*[(freq.mean(), freq.std()) for freq in np.array(lista_d50).transpose()])
#axs[-1].plot(freqs*1000, media_d50, ls='-')
#axs[-1].plot(freqs*1000, np.array(lista_d50).transpose()) # plota todas as respostas em um só gráfico
axs[-1].errorbar(freqs*1000, media_d50, yerr=desvio_d50, label='D50 com desvio padrão')
axs[-1].set_title('D50 médio da igreja')
axs[-1].fill_between(freqs*1000, .3, .7, color='green', alpha=.2, label='Faixa de variação típica (0.3 à 0.7)')
# escala logarítimica
axs[-1].set_xscale('log', basex=2)
axs[-1].set_xticks(1000 * freqs)
axs[-1].xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
# linha mediana
mediana = np.median(media_d50)
d50_1k = media_d50[12]
axs[-1].axhline(y=d50_1k, color='red', ls='--', label='D50 em 1.000 Hz = {:.2f}'.format(d50_1k))
axs[-1].axis([min(freqs) * 1000, fs*1000/2, 0, 1]) # limites
# legenda
handles, labels = axs[-1].get_legend_handles_labels()
axs[-1].legend(handles, labels,loc='lower right')
# insere o texto dos eixos
fig.text(.5, .001, 'Frequency (Hz)', ha='center')
fig.text(-.02, .5, 'D50', va='center', rotation='vertical')
#fig.savefig('D50.pdf', format='pdf')
O centróide espectral é uma indicação do "centro de massa" de um sinal. É uma média das frequências do sinal, ponderada por suas amplitudes, i.e., $$\text{Centroide} = \frac{\sum_{n=0}^{N-1}f(n)X(n)}{\sum_{n=0}^{N-1}X(n)}$$ onde f(n) é a n-ésima frequência central dada pela FFT e X(n) sua amplitude.
def centroid(ir, fs):
X = np.abs(rfft(ir)) # magnitudes
#X = X / sum(X) # normalização
freqs = np.linspace(0, 1, len(ir))
bins = np.abs(np.fft.fftfreq(len(ir), 1.0 / fs)) # frequencias da fft
return np.sum(bins * X) / np.sum(X) #* bins[-1] # média ponderada
# cria uma lista das bandas padrão de 1/3 de oitava em kHz
freqs = np.array([.020, .025, .032, .040, .050, .063, .080, .100, .125, .160, .200, .250, .315, .400, .500, .630, .800,
1.0, 1.25, 1.6, 2.0, 2.5, 3.15, 4.0, 5.0, 6.3, 8.0, 10.0, 12.5, 16.0])
# cria gráficos suficientes para plotar todas as respostas
fig, graficos = plt.subplots(num_respostas + 1, figsize=(15,20)) # 'graficos' é uma lista de gráficos
plt.tight_layout() # apenas estético
lista_centroides = []
# plota cada resposta nos graficos seguintes
for indice, resposta in enumerate(db_fft_respostas): # para cada resposta em frequência
graficos[indice].plot(f, resposta) # plota o gráfico da resposta
graficos[indice].set_title('Resposta {}'.format(indice + 1)) # título da resposta
graficos[indice].axis([16, (fs*1000)/2, 0, 60]) # ajusta os limites dos eixos
# escala logarítmica
axs[indice].set_xscale('log', basex=2)
axs[indice].set_xticks(1000 * freqs)
axs[indice].xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
# centróide espectral
centroide = centroid(resposta, fs*1000)
lista_centroides.append(centroide)
graficos[indice].axvline(centroide, color='red', ls='-', label='Centróide: {:.0f} Hz'.format(centroide))
# legenda
handles, labels = graficos[indice].get_legend_handles_labels()
graficos[indice].legend(handles, labels)
# plota a média dos centróides
centroide_medio, desvio_centroide = np.mean(lista_centroides), np.std(lista_centroides)
resposta_media = sum(db_fft_respostas) / len(db_fft_respostas)
graficos[-1].plot(f, resposta_media)
graficos[-1].set_title('Centróide de frequência médio da igreja')
graficos[-1].axis([16, (fs*1000) / 2, 0, 60])
graficos[-1].axvline(centroide_medio, color='red', ls='-', label='Centróide: {:.0f} Hz'.format(centroide_medio))
handles, labels = graficos[-1].get_legend_handles_labels()
graficos[-1].legend(handles, labels)
# insere o texto dos eixos
fig.text(.5, .001, 'Frequency (Hz)', ha='center')
fig.text(-.02, .5, 'Amplitude (dB)', va='center', rotation='vertical')
#fig.savefig('centroide_de_frequencia.png', format='png')
rate, voice = audioread('Voice.16bit.wav')
_, musica_lenta = audioread('HarpMono.wav')
_, musica_rapida = audioread('musica_rapida.wav')
auralizacao = np.convolve(voice, irs[1][:int(2e3*fs)])
titulo1 = IPython.display.Markdown(data='Original') # título
audio1 = IPython.display.Audio(data=voice, rate=rate) # áudio
IPython.display.display(titulo1)
IPython.display.display(audio1)
titulo = IPython.display.Markdown(data='Auralização') # título
audio = IPython.display.Audio(data=auralizacao, rate=rate) # áudio
IPython.display.display(titulo)
IPython.display.display(audio)
audiowrite(filename='auralizacao.wav',
rate=rate,
data=(0.9*auralizacao).astype(np.int16))
auralizacao_musica_lenta = np.convolve(musica_lenta, irs[1][:int(2e3*fs)])
auralizacao_musica_rapida = np.convolve(musica_rapida, irs[1][:int(2e3*fs)])
titulo1 = IPython.display.Markdown(data='Música 1 - Original') # título
audio1 = IPython.display.Audio(data=musica_lenta, rate=rate) # áudio
IPython.display.display(titulo1)
IPython.display.display(audio1)
titulo2 = IPython.display.Markdown(data='Música 2 - Original') # título
audio2 = IPython.display.Audio(data=musica_rapida, rate=rate) # áudio
IPython.display.display(titulo2)
IPython.display.display(audio2)
titulo3 = IPython.display.Markdown(data='Música 1 - Auralização') # título
audio3 = IPython.display.Audio(data=auralizacao_musica_lenta, rate=rate) # áudio
IPython.display.display(titulo3)
IPython.display.display(audio3)
audiowrite(filename='auralizacao_musica_lenta.wav',
rate=rate,
data=(0.9*auralizacao_musica_lenta).astype(np.int16))
titulo4 = IPython.display.Markdown(data='Música 2 - Auralização') # título
audio4 = IPython.display.Audio(data=auralizacao_musica_rapida, rate=rate) # áudio
IPython.display.display(titulo4)
IPython.display.display(audio4)
audiowrite(filename='auralizacao_musica_rapida.wav',
rate=rate,
data=(0.9*auralizacao_musica_rapida).astype(np.int16))